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Abstract 

First passage distributions of semi-Markov processes are of interest in fields such as 
reliability, survival analysis, and many others. The problem of finding or computing 
first passage distributions is, in general, quite challenging. We take the approach of us- 
ing characteristic functions (or Fourier transforms) and inverting them, to numerically 
calculate the first passage distribution. Numerical inversion of characteristic functions 
can be numerically unstable for a general probability measure, however, we show for 
lattice distributions they can be quickly calculated using the inverse discrete Fourier 
transform. Using the fast Fourier transform algorithm these computations can be ex- 
tremely fast. In addition to the speed of this approach, we are able to prove a few 
useful bounds for the numerical inversion error of the characteristic functions. These 
error bounds rely on the existence of a first or second moment of the distribution, or on 
an eventual monotonicity condition. We demonstrate these techniques in an example 
and include R-code. 

KEY WORDS: characteristic function, discrete Weibull, first passage dis- 
tribution, FAST Fourier transform, semi-Markov process, statistical flow- 
graph 

1 Introduction 

Statistical literature abounds with proofs using characteristic functions (CFs). Asymptotic 
results such as the central limit theorem rely heavily on the properties of CFs. However, 
when it comes to applied statistics the reverse is true. In general statisticians seem very 
uncomfortable using the CF or numeric approximations of it, even though they routinely 
use numerical approximations and calculations for various other procedures. This could 
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be partly due to the fact that CFs are complex functions, but this is a deterrent based on 
fear of the unknown, not difficulty. 

This is not the case in many applied sciences. Engineers make heavy use of the Fourier 
transform and are quite comfortable using it in practice. The Fourier transform and the 
characteristic function are the same function with a change of variables. To use one over 
the other is a matter of preference, the curves, as drawn, are identical. There are some 
areas in applied statistics where use of CFs are helpful. One such area is semi-Markov 
processes (SMPs), specifically find ing first pas sage distributions from one state to another, 
as in statistical flowgraph models iHuzurbaza r (2005). Using characteristic functions is a 
natural way to find first passage distributions in SMPs. 

For example consider the graphical representation of a SMP in FigurelU where the nodes 
on the graph represent a state, and the branches represent waiting time distributions before 
transition to the next state. Each branch is labeled with a characteristic function (f^j^s) 
and a probability of taking that branch. The characteristic function of the first passage 
distribution is: 

(1 -P)¥'12(S)V'23('5) 



(1) 



the 



1 - Pfl2 is)ip21 (S) 

denotes a first passage. For details on finding the first pas sage CF see Huzurbazail 



(|2005l ) or for a more theoretical development refer to iPvkd (|l96lh . Once the CFs for the 
branches are known, it is trivial to calculate the CF of the first passage from state 1 to 
3. The challenge is to convert the CF to a distribution. In most cases this must be done 
numerically. 
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Figure 1: An example semi-Markov process. 



Work has been done to show that the discrete Fourier transform (DFT) can accurately 
invert characteristic functi ons given suffic iently fast decay of both the density and the 
characteristic function (see Hughett ( 19981 )). One class of distributions that has not been 
well studied for inversion using the DFT is discrete distributions. Specifically we focus 
on distributions with support nAt for n = 0,=bl,±2, ... and At € M"*". This type of 
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distribution has been termed a lattice or an arithmetic distribution, these include most 
of the well studied discrete distributions such as the binomial, Poisson, negative binomial 
and many others. 

The primary goal of this paper is to promote the use of the inverse discrete Fourier 
transform (iDFT) for the inversion of CFs to probability mass functions (PMFs). A sec- 
ondary benefit of using the DFT is easy calculation of CFs when they do not have a closed 
form. 

This paper is outlined as follows. The next section details some properties of the DFT, 
in particular how the DFT is applied when dealing with lattice distributions. Next, the 
DF T is shown in action, in an applica tion. We include code using the R software package 
(see R Development Core Team ( 20091 )). Finally, we conclude with a few comments. 



2 The discrete Fourier transform 

For a random variable T, with distribution F, its characteristic function is 

fOO 

e*^* dF. (2) 



oo 



Let /(w) be the Fourier transform of T, then the relationship between the characteristic 
function and the Fourier transform is 

ip^{-27Tu;) = f{uj). (3) 

In the reminder of the paper we use Fourier transforms, with the understanding that the 
CF is easily obtained from the FT. We do this primarily for ease in applying DFT methods, 
specifically using the fast Fourier transform (FFT) algorithm. 

The DFT can be used to approximate the Fourier transform and its inverse. Under 
certain conditi ons, the error of this approx imation can be controlled. An excellent reference 



on the DFT is iBriggs and HensonI (jl995l ). To use the DFT, the support of the RV must 
be discretized at evenly spaced points, where At is the space between adjacent points. We 
term these points the support of the DFT samples. Similarly the support of the iDFT 
samples the FT at evenly spaced intervals of width Alo in the transform domain. These 
sampling widths are related by the reciprocity relation 

AtAu = ^. (4) 
With this information we apply the DFT and iDFT to lattice distributions. 
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2.1 DFT of lattice distributions 



In this section error bounds for lattice distributions with nonnegative support are derived. 
If f{tn) is the PMF of a lattice distribution, then the DFT of that distribution is: 



N-l 

Fk= ^/(tn)e-'"'^, for fc = 0,l,...,7V-l. (5) 

n=0 

Clearly the utility of the DFT as an approximation of the FT is dependent on the size of A'^. 
The major benefit of using the DFT with discrete random variables is that unlike RVs with 
PDFs, when sampling from lattice distributions we are able to capture essentially all the 
information about the random variable. In this paper we use the word "sampling" loosely 
to mean the PMF values of the support that are included in the DFT. When sampling a RV 
with a PDF the measure of our finite number of samples is zero, whereas when sampling 
from a certain class of lattice distributions our "sample" has measure that can approach 
one. In essence, we are losing little or no information when sampling from a DRV, or in 
other words, with a finite number of points we obtain more of a "census" than a "sample" 
of where the probability is located in a distribution. This allows the DFT to accurately 
approximate the FT at specified points given N sufficiently large. 

We confine our focus to nonnegative lattice RVs T where, given any e > 0, there exists 
an such that X^^q^ P{T = ti) > 1 — e. Nearly all properly defined nonnegative lattice 
RVs fall into this class. There may be a few pathological cases which do not, such as the 
counting measure on all the integers, however, in many cases one would argue that these are 
not properly defined distributions. The assumption that T is nonnegative is only slightly 
less general than the unrestricted case, and derivations for lattice RVs with unbounded 
support would be very similar. 

To find the error bound for the DFT on lattice distributions one only needs to look at 
the definition of the Fourier transform. Let /(in) = -P(T = n At), then the definition of 
the Fourier transform is: 

oo oo 

f{Lj) = J2 /(i„)e-^™""^* = Y^P{T = nAt)e-2™"^* (6) 

71=0 n=0 

Now our goal is to bound the error |/(wfc) — F^l for all index k = 0,1, . . . , N — 1. 
Therefore we must find an optimal Uk that matches well with the DFT's F^. To do this 
we simply set the exponents of Equations [5] and [6] equal to each other and solve for lo^ ■ 
Therefore 

nk 

— 27riuJi^nAt = — 27ri — , (7) 

implies that 

wfc = ir^- (8) 
AtN ^ ' 
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Now using the constraint that 



N-l 



N-1 



n=0 n=0 

we show that the error can be controlled. 

Lemma 2.1 For a nonnegative lattice random variable, defined on the support nAt for 
n = 0, 1, . . . , oo, the pointwise error of the forward DFT is less than or equal to P{T > 
NAt). 



Proof 



f{uk) - Fk 



N-l 



P{T = nAt)e-'^''^^''^* -Y^P{T = nAt)e 

n=0 n=0 
oo AT-l 

P{T = nAt)e-2'^^t _ ^ p(r = nAt)e-^''''^ 

n=0 

oo 

Y P{T = nAt)e 



n=0 



n=N 

oo 



e N 



n=N 

oo 



£p(r = nAt) 



n=N 



J2fitn) =P{T>NAt). I 



n=N 

To give an example, consider a Poisson random variable U, with rate parameter A = 5. 

The FT of U is f{uj) = exp{5e-2'^*'^ - 5}. If we choose N = 16, where At = 1, then 
UJk = k/16 for A: = 0, 1, ... , 15. Wc set the error e = P{U > 15) ^ 0.000069. The greatest 
error occurs when computing |/(0/16) — Fo| « |1 — 0.999931| = 0.000069. Therefore we 

can see in this example that f{ojk) — Fk < s. 

Now that we have shown that the DFT can accurately approximate the Fourier trans- 
form of a DRV, at specified points, we focus on the inverse DFT. This is the more difficult 
of the two problems, but is in general more useful. Given the Fourier transform of a DRV 
we now demonstrate how to accurately calculate the PMF. 

Given that f{uj) is a known function, we show that the inverse DFT has the same error 
bound as the forward DFT. The inverse DFT is defined as: 



(9) 



n=0 
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Lemma 2.2 Given the Fourier transform of a nonnegative lattice RV, T, evaluated at 
points k/{N At ), for k = 0, 1,...A^ — 1, the pointwise error of the inverse DFT is less than 
or equal to P{T >N At). 

We introduce a couple of notational items to help with the following proof. Let "mod" be 
the modulo operator, so for the positive integer b, and non-negative integers a, n and r; 
then r = a mod b, where n is the largest non-negative integer that satisfies the equation 
a = r + nb. This implies that r < b. Also, let 

_ J 1 if "condition" is true , . 

J(condition) - I Q otherwise 

We'll call /(condition) f^e indicator function. 

Proof Using Equations and [S] we can write the iDFT as: 

JV-l oo 



n=0 j=0 
N-l oo 

n=0 j=0 

oo N-l 

j=0 n=0 
^ oo 

= nY.^^^ = ■?'^*) ^ kk=i mod N) 



j=0 

N-l oo 

= P{T = J At) j^od JV) + E = -^^^^ ^('^=i N) 

j=0 j=N 

oo 

= P{T = kAt) + ^ P(T = J At) ^od AT). 
Therefore we can now express the error as: 

oo 

\fk - P{T = kAt)\ = P{T = kAt) + Y,P{T = j/^t)\k=jmodN)-P{T = kAt) 

j=N 

< P{T > NAt). I 

This result is not surprising, because of the close relationship between the forward and 
inverse DFT. 
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The next result assumes the the FT is unknown and must be approximated by the 
DFT, following which, the iDFT is used to invert the approximated FT to a PMF. This 
theorem shows this error is also bounded. 



Theorem 2.3 Given that the FT of a nonnegative lattice random variable is known up 
to a fixed error bound, then the error bound of the iDFT of this FT can be bounded by 
2P{T >NAt). 

Proof Mimicking the proof of Lemma 12.21 we replace f{oJn) with f{oJn) + where e = 
P{T > NAt) and |(5„| < e. This gives: 



fk 



N-l 



n=0 



2lTiu}r, kAt 



N-l 



2ni!- 



n=0 n=0 



= P{T = kAt) + Y,PiT = j^t)l^k=j mod N) 

Therefore if we only have an approximation for f{u)n) then, 



N-l 



\fk - PiT = kAt)\ 



n=0 



oo N-l 
J2 PiT = jAt) \k=j mod N) + J^Y. 



j=N 

oo 



X] = jAt) l(f,=j mod N) 

j=N 

n=0 

N-l 
£ \ - 

< e + —yi 



+ 



n=0 
N-l 



1 \ ^ r- l-jri — 



n=0 



< e 



e N 



2e. 



n=0 



What we have not addressed, thus far, is the error introduced when manipulating FTs 
in the Fourier domain. The primary reason we use CFs and FTs is for the convenient 
theoretical properties but when f{u>n) is estimated, then when mathematical manipulation 
such as multiplication or division occurs the error is magnified. We do not specifically 
address this problem, but warn the practitioner of this added source of error. Additionally, 
machine precision error also plays a part in the overall error bound. 
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2.2 Finding the inversion error bounds 

The primary problem when using the inverse DFT is that one may not know a priori N 
such that 

Af-l 

J2 P(.T = nAt) >l-e, 

n=0 

for a given e. So although the error is bounded, this bound is not obvious without the 
exact PMF, which we do not usually have. There are a few approaches to overcome this 
difficulty. One solution is to use an inequality to bound P(T > nAt). One bound is 
Markov's inequality: 

P{T >a)< for a > 0. 

a 

Therefore if the first moment exists for the RV T, and we either know or can estimate it, 
this method provides a conservative error bound. The primary issue is finding E{T). In our 
case the random variable T's distribution is a first passage distribution in a semi-Markov 
process. It turns out if all the expectations to the individual transitions in the SMP are 
known, then finding the first passage expectation is the solution to a set of linear equations. 
Literature exist s that explains h o w one can find the mo ments of first passage distributions, 
for example see IZhang and Houl (I2OI2I ') and|^ (jlQSSl ). 



Another inequality that can assist in bounding the error is Cantelli's inequality. For a 
random variable T with ji = E{T), = Var(T), and A; > we have 

P{T >ka + fi)< ^ 



1 + A:2- 



Therefore if Var(T) < 00, and given E{T) and Var{T), Cantelli's inequality can be useful 
in bounding the DFT error. Where k > cr/fi, Cantelli's inequality provides a smaller 
bound and should be used over Markov's inequality. Finding the va riance of a first passage 
random variable is similar to the expectation, again see lYa^ (| 1985 11 for details. 



If the two inequalities prove to be impractical, then another approach can be taken if 
certain assumptions can be made about the distribution of T. 

Theorem 2.4 If for some index M where the following inequalities are satisfied P{T = 
nAt) < P{T = mAt) and M < m < n, then the error for the iDFT can be expressed 
as: Choose an even N (the number of DFT samples) such that N > 2M, then for all 
n = 0,l,...,{N/2- 1) 

|/„-P(T = nAt)j </^/2+„. 

Essentially what this is saying is that after a certain point, if the first passage PMF, f{t), 
is monotonically nonincreasing, then we can bound the error using that fact alone. This 
technique is certainly not universally applicable, but it does provide a working error bound 
if monotonicity eventually occurs. To our knowledge PMFs used for modeling processes 
will often meet this assumption. 
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Proof From the proof of Lemma 12.21 we know 

oo 

fn = ^ P{T = jAt) I(n=j mod N)i 
3=0 

Which can be rewritten as 

oo 

fn = Y.P{T = {n + jN)At). 

j=0 

Which imphes 

oo 

fN/2+n = J2p{T= {N/2 + n + jW) At) 

therefore 

oo 

|/„ - P{T = nAt) I = ^ P(r = (n + jN)At) 

but by monotonicity we know 

P{T= {N/2 + n + jN) At) >P{T = {n+{j + l)N) At) for all j = 0, 1, 2 . . . 
This implies 

oo 

Therefore 

|/„-P(r = nAt)| </^/2+„. I 

This bound can be the more effective then the previous bounds, if the condition exists 
for it to be applied. However, it is difficult to prove the monotonicity condition in a 
particular situation, and if sufficient doubt exists, another error bounding method should 
be applied. 

Before we move to the application section, we suggest an even broader condition than 
eventual monotonicity that will provide the same error bound. This condition is: 

Theorem 2.5 If there exists some nonnegative integers n and N/2 and a positive integer 
k, such that, ifn < N/2, and P{T ={n + kN/2) At) >P{T={n + {k + l)N/2) At) then, 
for a DFT with N samples, 

|/„-P(r = nAt)| </„+^/2. 

These types of distributions might be termed eventually periodic nonincreasing. This con- 
dition is very broad and can be applied in nearly all first passage distributions, although 
we are sure there exist some pathological cases that do not fit in this class of distributions. 
The proof of Theorem 12.51 follows the same argument as Theorem 12.41 

In the next section we demonstrate using the DFT to calculate FTs, their inversion, 
and also it's error bounds. 
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3 An application of discrete time semi-Markov processes 



Consider the problem described in Barbu and Liminios ( 20081 ). the waste treatment for a 
textile factory. If the treatment facility is working then the waste can be discarded; if the 
unit fails the untreated waste must be stored in a holding tank until the treatment facility 
is again functioning. However, if the repair to the treatment facility takes too long the 
holding tank becomes full, and the textile factory must halt production until the repair is 
complete. Figure [2] is a graphical depiction of the process. 




Working 



Cp3i(s) 



912(5) 





(1-P)CP23(S) 



Treatment failed, 
tank not full 



Treatment failed, 
tank full 



Figure 2: A graphical depiction of the textile waste treatment process. 



We use the model parameterization of iBarbu and LiminiosI (|2008l ). which is: 



T12 ~ geometric(0.8), 

T21 ~ discrete Weibull(0. 3, 0.5), 

T23 ~ discrete Weibull(0.5, 0.7), and 

T31 ~ discrete Weibull(0. 6, 0.9). 



The discrete Weibull as defined in lNakagawa and Osakil (jl975l ) (slightly modified) is 



f{t\q,b) 



f iftG{l,2,...} 
otherwise. 



The SMP parameter p (see Figure [2]) is set at 0.95, or in words, on average 95% of the 
repairs to the treatment facility are completed before the tank becomes full. 

A few of the quantities that may be of interest from this process are the first passage 
PMF from state 1 to state 3, the expected amount of time spent in state 3 at some time t, 
and the probability the factory is shutdown due t o waste treatment failure at some time t. 
The formulas for these quantities can be found in lWarr and Collins For this paper 

we focus on finding the first passage PMF from state 1 to state 3. 
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Table 1: The mean and variance of waiting time random variables in the textile waste 
treatment process. 





Mean 


Variance 




1.2500 


0.3125 


T21 


2.0769 


9.0765 


T23 


2.7300 


9.4618 




67.191 


4394.1 



The first passage CF is given in Equation [H where we replace CFs with FTs in all 
instances. So to calculate the first passage FT from state 1 to state 3, we need the FT of 
T12, T21, and T23. The FT of geometric RV T12 is known in closed form, however, the FT 
for the discrete Weibull is not. So in our computations we use the exact FT for T12 and 
use the DFT approximation for the other two. 

We first want to calculate the mean and variance of T12 , T21 , and T23 which will allow 
us to use the Cantelli inequality to bound our error. Table [T] shows the expectations and 
variances for these RVs. These quantities allow us to calculate the number we need 
to use to obtain a desired error bound. We choose e = 10~^ which implies k > 66356. 
Therefore if A'^ = 2^^ the pointwise error for each probability will be smaller than 2e. 
The computation time of calculating the first passage PMF is less than one second (on a 
Windows 7 PC with an AMD Athlon'^'^ 7750 2.7 GHz dual -core processor). 

The plot in Figure [3] shows what the first passage PMF from state 1 to state 3 looks like 
for the first one hundred points. Using the developed theory in this paper, we are assured 
the pointwise error is less than or equal to 2e. Again this ignores the error introduced by 
multiplying and dividing the FTs (as in Equation [T]) and the machine round-off error. 

Figure [3] indicates that the PMF of this first passage distribution is monotonically 
decreasing after roughly time 20. Therefore we can attempt to apply Theorem 12. 4[ Since 
we estimated that P(Ti3* = 609) < e Therefore if we choose A^ = 1218 we only need to 
calculate roughly 1% of the values in the FFT and can still claim the pointwise error is 
less than or equal to 2e. This example shows the potential computation savings of using 
Theorem 12.41 for first passage random variables. Also, Theorem 12.41 does not assume a 
finite mean or variance, but a monotonically decreasing PMF after a selected point on the 
support. In this example we worked our way, by trial and error, from a large A^ down to 
a smaller A^, however, we are also able to work our way up from a small A^ to a larger A^. 
This can reduce the computation time to find the desired error bound. 
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Figure 3: The PMF of the first passage from state 1, "Working", to state 3, "Treatment 
failed, tank full". 



4 Discussion 

It has been our goal to show for lattice random variables the DFT/FFT computes the 
desired first passage PMF extremely quickly and to a prescribed accuracy. It is often very 
desirable to handle probability problems in the transformed CF domain, but computational 
methods to invert them back into the time or spatial domain have been the largest deterrent. 
The DFT is a excellent solution for the inversion of characteristic functions to PMFs. The 
DFT is very fast using the FFT algorithm and it also provides error bounds for estimates. 
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6 Appendix — R Code 



##Discrete Weibull PMF 
dweibulldisc<-f unction (x,q,b){ 

templ=q"((x-l)"b)-q"(x"b) ; temp2=x-f loor (x) 

tenip3=tenip2<le-15 ; temp4<-x==0; tempi [temp4] <-0 

return (templ*temp3) 

} 

##p is the probability of the treatment facility will be fixed 
## before the holding tank becomes full 
p = 0.95 

N = 2'17 

support = 0: (N-1) 

## The DFT samples of the two of the transition distributions 
T21 = dweibulldisc (support , . 3 , . 5) 
T23 = dweibulldisc (support, 0.5, 0.7) 

## The Fourier transforms of the 3 transition distributions 

F12 = p*exp(-2*pi*li*support/N)/(l-(l-p)*exp(-2*pi*li*support/N)) 

F21 = fft(T21); F23 = fft(T23) 

## The Fourier transform of the first passage distribution from 

## state 1 to state 3 

F_lst_Pas_13 = ((l-p)*F12*F23)/(l-p*F12*F21) 

## The first passage PMF from state 1 to state 3 
T_lst_Pas_13 = l/N*Re( f f t (F_lst_Pas_13, inverse=T) ) 
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